#____________________________________________________________#
### Paper: Sex Crime and Punishment in the #MeToo Era 
### Authors: Susanne Schwarz, Dara Kay Cohen, Matthew Baum 
### Code: Replication Code for Analyses and Graphs in Paper 
#____________________________________________________________#

rm(list = ls())

library(foreign)
library(multcomp)
library(ggplot2)
library(cjoint)
library(plyr)
library(dplyr)
require(readr)
library(commarobust)
library(readstata13)
library(knitr)
library(kableExtra)
library(tidyverse)
library(estimatr)
library(stringr)
library(irr)

#____________________#
# 1) SET-UP 
#____________________#

#Load Data

load("~/main.study.data.RData")
load("~/follow.up1.clean.RData")
load("~/coded.data.follow.up1.RData")



### Create Subsets for Different conditions 

d.rep.rap <- subset(d, d$RAPE==1 & d$REPORT==1)       #Rape Reporting
d.rep.rob <- subset(d, d$ROBBERY==1 & d$REPORT==1)    #Robbery Reporting
d.pun.rap <- subset(d, d$RAPE==1 & d$PUNISH==1)       #Rape Punishment 
d.pun.rob <- subset(d, d$ROBBERY==1 & d$PUNISH==1)    #Robbery Punishment
d.report <- subset(d, d$REPORT==1)                    #Reporting Overall
d.punish <- subset(d, d$PUNISH==1)                    #Punishment Overall 



# Plot Prep
LABELS <- c("Victim: Female (vs. Male)",
            "Victim: Black (vs. White)",
            "Victim: Felony Record (vs. No Record)",
            "Victim: Multiple Partners (vs. Married)",
            "Victim: Single (vs. Married)",
            "Perpetrator: Black (vs. White)",
            "Perpetrator: Athlete (vs. Worker)",
            "Perpetrator: Business man (vs. Worker)",
            "Perpetrator: Acquaintance (vs. Stranger)",
            "Location: Park (vs. Home)", 
            "Location: Party (vs. Home)",
            "Attire: Club Outfit (vs. Work Outfit)")


#import helper functions
source("~/helper.functions.R")


#____________________________________________________________#
# 2) ANALYSIS: AMCE Estimation using OLS regression
#____________________________________________________________#


# a) Rape Reporting

lm.model1 <- lm_robust(CHOSEN ~ VIC_SEX + VIC_ETHN + VIC_CRIMINAL + VIC_MARITAL + PERP_ETHN + PERP_PROF + PERP_RELATION + LOCATION + OUTFIT, 
                data=d.rep.rap, 
                clusters=ID,
                weights=WEIGHT) #AMCE estimation; SE clustered on respondent level


lm.model1.cluster <- tidy(lm.model1)  


# b) Rape Punishment

lm.model2 <- lm_robust(CHOSEN ~ VIC_SEX + VIC_ETHN + VIC_CRIMINAL + VIC_MARITAL + PERP_ETHN + PERP_PROF + PERP_RELATION + LOCATION + OUTFIT, 
                data=d.pun.rap, 
                cluster=ID,
                weights=WEIGHT) #AMCE estimation; SE clustered on respondent level

lm.model2.cluster <- tidy(lm.model2) 


# c) Figure 1: Main Paper 

model1.plot<- lm.model1.cluster[-1,]
model2.plot<- lm.model2.cluster[-1,]


# prepare output for ggplot
model1.plot <- model1.plot  %>% mutate(model = "Reporting",
                                       term = LABELS)
model2.plot <- model2.plot %>%mutate(model = "Punishment",
                                     term = LABELS)

gg_df <- bind_rows(model1.plot, model2.plot) %>%
  mutate(model = factor(model, levels = c("Reporting", "Punishment")))

# plot

ggplot(data = gg_df, aes(estimate, term)) +
  geom_point() +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0, linetype=ifelse(gg_df$conf.low>0 | gg_df$conf.high<0, "solid", "dashed")) +
  facet_wrap( ~ model) +
  geom_vline(xintercept = 0) +
  xlab('Change in Pr(Case Chosen)') + 
  ylab('Profile Attributes') +
  theme_bw() +
  theme(axis.title.y = element_blank(),
        axis.title.x = element_text(size=14),
        axis.text.y = element_text(size = 14),
        axis.text.x = element_text(size = 14),
        strip.text = element_text(size=14)) +
  ggtitle("Rape Reporting & Punishment") 


#____________________________________________________________#
# 3) ANALYSIS: Predicted Probabilities of Being Chosen
#____________________________________________________________#

#a) Rape reporting 

preds <- as.data.frame(predict(lm.model1, newdata=d.rep.rap, se.fit = TRUE))

d.rep.rap <- cbind(d.rep.rap, preds)
d.rep.rap <- d.rep.rap %>% 
  mutate(PCT = ntile(fit, 100),
         QUT = ntile(fit,5))

quantiles <- list(1, 25, 50, 75, 99)

x <- subset(d.rep.rap, d.rep.rap$PCT %in% quantiles, select=c("VIC_SEX", "VIC_ETHN", "VIC_CRIMINAL", "VIC_MARITAL", "PERP_ETHN",
                                                              "PERP_PROF", "PERP_RELATION","LOCATION", "OUTFIT", 
                                                              "CHOSEN", "fit", "se.fit", "PCT", "QUT"))

new_df <- x %>% group_by(PCT) %>% sample_n(1)

new_df$term <- NA
for(i in 1:nrow(new_df)){
  new_df$term[i] <- getlabels(new_df[i,colnames(new_df)[1:9]])
}

# Add Confidence Intervals 

new_df$ui <- new_df$fit+1.96*new_df$se.fit
new_df$li <- new_df$fit-1.96*new_df$se.fit


# Figure 2: Panel A

ggplot(new_df, aes(fit, y= reorder(term, -fit)), label = c("1st Percentile", "25th Percentile", "50th Percentile", "75th Percentile", "99th Percentile")) +
  geom_point() +
  geom_text(aes(label = c("1st Percentile", "25th Percentile", "50th Percentile", "75th Percentile", "99th Percentile")), nudge_y = 0.2) +
  geom_errorbarh(aes(xmin = li, xmax = ui), height = 0, width=1) +
  scale_y_discrete(labels = function(x) str_wrap(x, width = 60)) +
  geom_vline(xintercept = 0.5) + 
  xlab('Pr(Case Chosen for Reporting)') + 
  ylab('PCase Profiles') +
  theme_bw() +
  theme(axis.title.y = element_blank(),
        axis.title.x = element_text(size=11),
        axis.text.y = element_text(size = 14),
        axis.text.x = element_text(size = 11),
        strip.text = element_text(size=14),
        plot.title = element_text(size = 16, face = "bold")) +
  ggtitle("Panel A: Reporting") 



# b) Rape Punishment 

preds <- as.data.frame(predict(lm.model2, newdata=d.pun.rap, se.fit = TRUE))

d.pun.rap <- cbind(d.pun.rap, preds)
d.pun.rap <- d.pun.rap %>% 
  mutate(PCT = ntile(fit, 100),
         QUT = ntile(fit,5))

quantiles <- list(1, 25, 50, 75, 99)

x <- subset(d.pun.rap, d.pun.rap$PCT %in% quantiles, select=c("VIC_SEX", "VIC_ETHN", "VIC_CRIMINAL", "VIC_MARITAL", "PERP_ETHN",
                                                              "PERP_PROF", "PERP_RELATION","LOCATION", "OUTFIT", 
                                                              "CHOSEN", "fit", "se.fit", "PCT", "QUT"))

new_df2 <- x %>% group_by(PCT) %>% sample_n(1)

new_df2$term <- NA
for(i in 1:nrow(new_df2)){
  new_df2$term[i] <- getlabels(new_df2[i,colnames(new_df2)[1:9]])
}

new_df2$ui <- new_df2$fit+1.96*new_df2$se.fit
new_df2$li <- new_df2$fit-1.96*new_df2$se.fit


ggplot(new_df2, aes(fit, y= reorder(term, -fit)), label = c("1st Percentile", "25th Percentile", "50th Percentile", "75th Percentile", "99th Percentile")) +
  geom_point() +
  geom_text(aes(label = c("1st Percentile", "25th Percentile", "50th Percentile", "75th Percentile", "99th Percentile")), nudge_y = 0.2) +
  geom_errorbarh(aes(xmin = li, xmax = ui), height = 0, width=1) +
  scale_y_discrete(labels = function(x) str_wrap(x, width = 60)) +
  geom_vline(xintercept = 0.5) + 
  xlab('Pr(Case Chosen for Punishment)') + 
  ylab('PCase Profiles') +
  theme_bw() +
  theme(axis.title.y = element_blank(),
        axis.title.x = element_text(size=11),
        axis.text.y = element_text(size = 14),
        axis.text.x = element_text(size = 11),
        strip.text = element_text(size=14),
        plot.title = element_text(size = 16, face = "bold")) +
  ggtitle("Panel B: Punishment") 



#____________________________________________________________#
# 4) ANALYSIS: Rape vs. Robbery
#____________________________________________________________#


# a) Reporting condition

#RAPE
lm.modelA <- tidy(lm_robust(CHOSEN ~ VIC_SEX + VIC_ETHN + VIC_CRIMINAL + VIC_MARITAL + PERP_ETHN + PERP_PROF + PERP_RELATION + LOCATION + OUTFIT, 
                data=d.report[d.report$RAPE==1,], 
                cluster = ID,
                weights=WEIGHT))


#ROBBERY
lm.modelB <- tidy(lm_robust(CHOSEN ~ VIC_SEX + VIC_ETHN + VIC_CRIMINAL + VIC_MARITAL + PERP_ETHN + PERP_PROF + PERP_RELATION + LOCATION + OUTFIT,
                data=d.report[d.report$RAPE==0,], 
                cluster = ID,
                weights=WEIGHT))

# b) Punishment Condition 

# RAPE
lm.modelC <- tidy(lm_robust(CHOSEN ~ VIC_SEX + VIC_ETHN + VIC_CRIMINAL + VIC_MARITAL + PERP_ETHN + PERP_PROF + PERP_RELATION + LOCATION + OUTFIT, 
                data=d.punish[d.punish$RAPE==1,],
                cluster = ID,
                weights=WEIGHT))


# ROBBERY 

lm.modelD <- tidy(lm_robust(CHOSEN ~ VIC_SEX + VIC_ETHN + VIC_CRIMINAL + VIC_MARITAL + PERP_ETHN + PERP_PROF + PERP_RELATION + LOCATION + OUTFIT, 
                data=d.punish[d.punish$RAPE==0,], 
                cluster = ID,
                weights=WEIGHT))

# c) Figure 3: Main Paper 

# PANEL A
model.A.plot<- lm.modelA[-1,]
model.A.plot <- model.A.plot %>% mutate(
                           term = LABELS,
                           model = "Rape")

model.B.plot<- lm.modelB[-1,]
model.B.plot <- model.B.plot %>% mutate(
                          term = LABELS,
                          model = "Robbery")

allModelFrame.A <- rbind(model.A.plot, model.B.plot)

ggplot(data = allModelFrame.A, aes(estimate, term)) +
  geom_point() +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0, linetype=ifelse(allModelFrame.A$conf.low>0 | allModelFrame.A$conf.high<0, "solid", "dashed")) +
  facet_wrap( ~ model) +
  geom_vline(xintercept = 0) +
  xlab('Change in Pr(Case Chosen)') + 
  ylab('Profile Attributes') +
  theme_bw() +
  theme(axis.title.y = element_blank(),
        axis.title.x = element_text(size=14),
        axis.text.y = element_text(size = 14),
        axis.text.x = element_text(size = 14),
        strip.text = element_text(size=14), 
        plot.title = element_text(size = 16, face = "bold")) +
  ggtitle("Panel A: Reporting") 


# PANEL B

model.C.plot<- lm.modelC[-1,]
model.C.plot <- model.C.plot %>% mutate(
  term = LABELS,
  model = "Rape")

model.D.plot<- lm.modelD[-1,]
model.D.plot <- model.D.plot %>% mutate(
  term = LABELS,
  model = "Robbery")

allModelFrame.B <- rbind(model.C.plot, model.D.plot)

ggplot(data = allModelFrame.B, aes(estimate, term)) +
  geom_point() +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0, linetype=ifelse(allModelFrame.B$conf.low>0 | allModelFrame.B$conf.high<0, "solid", "dashed")) +
  facet_wrap( ~ model) +
  geom_vline(xintercept = 0) +
  xlab('Change in Pr(Case Chosen)') + 
  ylab('Profile Attributes') +
  theme_bw() +
  theme(axis.title.y = element_blank(),
        axis.title.x = element_text(size=14),
        axis.text.y = element_text(size = 14),
        axis.text.x = element_text(size = 14),
        strip.text = element_text(size=14), 
        plot.title = element_text(size = 16, face = "bold")) +
  ggtitle("Panel B: Punishment") 


# d) F-Tests 

#REPORTING

#Simple Model 
m1a <- lm(CHOSEN ~ VIC_SEX + VIC_ETHN + VIC_CRIMINAL + VIC_MARITAL + PERP_ETHN + PERP_PROF + PERP_RELATION + LOCATION + OUTFIT, 
          data=d.report, 
          weights=WEIGHT)

#Interaction model 
m1b <- lm(CHOSEN ~ VIC_SEX + VIC_ETHN + VIC_CRIMINAL + VIC_MARITAL + PERP_ETHN + PERP_PROF + PERP_RELATION + LOCATION + OUTFIT + RAPE + 
            RAPE*(VIC_SEX + VIC_ETHN + VIC_CRIMINAL + VIC_MARITAL + PERP_ETHN + PERP_PROF + PERP_RELATION + LOCATION + OUTFIT) , 
          data=d.report, 
          weights=WEIGHT)

#F-test of nested models
anova(m1a, m1b)

# PUNISHMENT 

#Simple Model 
m2a <- lm(CHOSEN ~ VIC_SEX + VIC_ETHN + VIC_CRIMINAL + VIC_MARITAL + PERP_ETHN + PERP_PROF + PERP_RELATION + LOCATION + OUTFIT, 
          data=d.punish, 
          weights=WEIGHT)

#Interaction model 
m2b <- lm(CHOSEN ~ VIC_SEX + VIC_ETHN + VIC_CRIMINAL + VIC_MARITAL + PERP_ETHN + PERP_PROF + PERP_RELATION + LOCATION + OUTFIT + RAPE + 
            RAPE*(VIC_SEX + VIC_ETHN + VIC_CRIMINAL + VIC_MARITAL + PERP_ETHN + PERP_PROF + PERP_RELATION + LOCATION + OUTFIT) , 
          data=d.punish, 
          weights=WEIGHT)
#F-test of nested models
anova(m2a, m2b)




#____________________________________________________________#
# 5) ANALYSIS: IE for Respondent Ideology
#____________________________________________________________#

# a) Rape Reporting

#Conservative 
lm.model.conservative.r <- tidy(lm_robust(CHOSEN ~ VIC_SEX + VIC_ETHN + VIC_CRIMINAL + VIC_MARITAL + PERP_ETHN + PERP_PROF + PERP_RELATION + LOCATION + OUTFIT, 
                                   data=d.rep.rap[d.rep.rap$CONSERV=="Conservative",], 
                                   clusters = ID,  
                                   weights=WEIGHT))


#Moderate
lm.model.moderate.r <- tidy(lm_robust(CHOSEN ~ VIC_SEX + VIC_ETHN + VIC_CRIMINAL + VIC_MARITAL + PERP_ETHN + PERP_PROF + PERP_RELATION + LOCATION + OUTFIT, 
                               data=d.rep.rap[d.rep.rap$CONSERV=="Moderate",], 
                               clusters= ID,
                               weights=WEIGHT))


#Liberal
lm.model.liberal.r <- tidy(lm_robust(CHOSEN ~ VIC_SEX + VIC_ETHN + VIC_CRIMINAL + VIC_MARITAL + PERP_ETHN + PERP_PROF + PERP_RELATION + LOCATION + OUTFIT, 
                              data=d.rep.rap[d.rep.rap$CONSERV=="Liberal",], 
                              clusters= ID,
                              weights=WEIGHT))


# b) Rape Punishment 

#Conservative
lm.model.conservative.p <- tidy(lm_robust(CHOSEN ~ VIC_SEX + VIC_ETHN + VIC_CRIMINAL + VIC_MARITAL + PERP_ETHN + PERP_PROF + PERP_RELATION + LOCATION + OUTFIT, 
                                   data=d.pun.rap[d.pun.rap$CONSERV=="Conservative",], 
                                   clusters = ID,  
                                   weights=WEIGHT))


#Moderate
lm.model.moderate.p <- tidy(lm_robust(CHOSEN ~ VIC_SEX + VIC_ETHN + VIC_CRIMINAL + VIC_MARITAL + PERP_ETHN + PERP_PROF + PERP_RELATION + LOCATION + OUTFIT, 
                               data=d.pun.rap[d.pun.rap$CONSERV=="Moderate",], 
                               clusters= ID,
                               weights=WEIGHT))


#Liberal
lm.model.liberal.p <- tidy(lm_robust(CHOSEN ~ VIC_SEX + VIC_ETHN + VIC_CRIMINAL + VIC_MARITAL + PERP_ETHN + PERP_PROF + PERP_RELATION + LOCATION + OUTFIT, 
                              data=d.pun.rap[d.pun.rap$CONSERV=="Liberal",], 
                              clusters= ID,
                              weights=WEIGHT))


# c) Plot Figure 4 in Main Paper

# PANEL A 
model.conserv.plot.r<- lm.model.conservative.r[-1,]
model.conserv.plot.r <- model.conserv.plot.r %>% mutate(
                          term= LABELS,
                          model = "01_Conservative")

model.moderate.plot.r<- lm.model.moderate.r[-1,]
model.moderate.plot.r <- model.moderate.plot.r %>% mutate(
                          term= LABELS,
                          model = "02_Moderate")

model.liberal.plot.r<- lm.model.liberal.r[-1,]
model.liberal.plot.r <- model.liberal.plot.r %>% mutate(
                          term= LABELS,
                          model = "03_Liberal")



allModelFrame <- rbind(model.conserv.plot.r, model.moderate.plot.r, model.liberal.plot.r)

ggplot(data = allModelFrame, aes(estimate, term)) +
  geom_point() +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0, linetype=ifelse(allModelFrame$conf.low>0 | allModelFrame$conf.high<0, "solid", "dashed")) +
  facet_wrap( ~ model) +
  geom_vline(xintercept = 0) +
  xlab('Change in Pr(Case Chosen)') + 
  ylab('Profile Attributes') +
  theme_bw() +
  theme(axis.title.y = element_blank(),
        axis.title.x = element_text(size=14),
        axis.text.y = element_text(size = 14),
        axis.text.x = element_text(size = 14),
        strip.text = element_text(size=14), 
        plot.title = element_text(size = 16, face = "bold")) +
  ggtitle("Panel A: Reporting") 


# PANEL B 
model.conserv.plot.p<- lm.model.conservative.p[-1,]
model.conserv.plot.p <- model.conserv.plot.p %>% mutate(
  term= LABELS,
  model = "01_Conservative")

model.moderate.plot.p<- lm.model.moderate.p[-1,]
model.moderate.plot.p <- model.moderate.plot.p %>% mutate(
  term= LABELS,
  model = "02_Moderate")

model.liberal.plot.p<- lm.model.liberal.p[-1,]
model.liberal.plot.p <- model.liberal.plot.p %>% mutate(
  term= LABELS,
  model = "03_Liberal")


allModelFrame <- rbind(model.conserv.plot.p, model.moderate.plot.p, model.liberal.plot.p)

ggplot(data = allModelFrame, aes(estimate, term)) +
  geom_point() +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0, linetype=ifelse(allModelFrame$conf.low>0 | allModelFrame$conf.high<0, "solid", "dashed")) +
  facet_wrap( ~ model) +
  geom_vline(xintercept = 0) +
  xlab('Change in Pr(Case Chosen)') + 
  ylab('Profile Attributes') +
  theme_bw() +
  theme(axis.title.y = element_blank(),
        axis.title.x = element_text(size=14),
        axis.text.y = element_text(size = 14),
        axis.text.x = element_text(size = 14),
        strip.text = element_text(size=14), 
        plot.title = element_text(size = 16, face = "bold")) +
  ggtitle("Panel B: Punishment") 


# d) F-Test of Nested Models 

#Note that the nested models below include additional respondent-level controls 

#Reporting 
modelA <- lm(CHOSEN ~ VIC_SEX + VIC_ETHN + VIC_CRIMINAL + VIC_MARITAL + PERP_ETHN + PERP_PROF + PERP_RELATION + LOCATION + OUTFIT + GENDER + AGE4 + RACETHNICITY + EDUC4 + REGION4, 
             data=d.rep.rap, 
             weights=WEIGHT)

modelB <- lm(CHOSEN ~ (VIC_SEX + VIC_ETHN + VIC_CRIMINAL + VIC_MARITAL + PERP_ETHN + PERP_PROF + PERP_RELATION + LOCATION + OUTFIT)*CONSERV + GENDER + AGE4 + RACETHNICITY + EDUC4 + REGION4, 
             data=d.rep.rap, 
             weights=WEIGHT)
anova(modelA, modelB)


#Punishment
modelC <- lm(CHOSEN ~ VIC_SEX + VIC_ETHN + VIC_CRIMINAL + VIC_MARITAL + PERP_ETHN + PERP_PROF + PERP_RELATION + LOCATION + OUTFIT + GENDER + AGE4 + RACETHNICITY + EDUC4 + REGION4, 
             data=d.pun.rap, 
             weights=WEIGHT)

modelD <- lm(CHOSEN ~ (VIC_SEX + VIC_ETHN + VIC_CRIMINAL + VIC_MARITAL + PERP_ETHN + PERP_PROF + PERP_RELATION + LOCATION + OUTFIT)*CONSERV + GENDER + AGE4 + RACETHNICITY + EDUC4 + REGION4, 
             data=d.pun.rap, 
             weights=WEIGHT)


anova(modelC, modelD)


#____________________________________________________________#
# 6) ANALYSIS: Main Results for Follow-Up #1
#____________________________________________________________#

d.report.fu1 <- subset(d.follow.up1, d.follow.up1$REPORT==1)
d.punish.fu1 <- subset(d.follow.up1, d.follow.up1$REPORT==0)

#Reporting
m1.fu <- tidy(lm_robust(CHOSEN ~ VIC_SEX + VIC_ETHN + VIC_CRIMINAL + VIC_MARITAL + PERP_ETHN + PERP_PROF + PERP_RELATION + LOCATION + OUTFIT, 
                       cluster=id,
                       data=d.report.fu1))

#Punishment 
m2.fu <- tidy(lm_robust(CHOSEN ~ VIC_SEX + VIC_ETHN + VIC_CRIMINAL + VIC_MARITAL + PERP_ETHN + PERP_PROF + PERP_RELATION + LOCATION + OUTFIT, 
                        cluster=id,
                        data=d.punish.fu1))


# Plot: Figure 5 in Main Paper 

model1.plot<- m1.fu[-1,]
model2.plot<- m2.fu[-1,]

model1.plot <- model1.plot  %>% mutate(model = "Reporting",
                                       term = LABELS)
model2.plot <- model2.plot %>%mutate(model = "Punishment",
                                     term=LABELS)

gg_df <- bind_rows(model1.plot, model2.plot) %>%
  mutate(model = factor(model, levels = c("Reporting", "Punishment")))

# plot
ggplot(data = gg_df, aes(estimate, term)) +
  geom_point() +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0, linetype=ifelse(gg_df$conf.low>0 | gg_df$conf.high<0, "solid", "dashed")) +
  facet_wrap( ~ model) +
  geom_vline(xintercept = 0) +
  xlab('Change in Pr(Case Chosen)') + 
  ylab('Profile Attributes') +
  theme_bw() +
  theme(axis.title.y = element_blank(),
        axis.title.x = element_text(size=14),
        axis.text.y = element_text(size = 14),
        axis.text.x = element_text(size = 14),
        strip.text = element_text(size=14)) +
  ggtitle("Rape Reporting & Punishment") 


#____________________________________________________________#
# 7) ANALYSIS: Inter Coder Reliability
#____________________________________________________________#


Coders <- sort(unique(coded.data$Coder))

data. <- lapply(1:length(Coders),function(i){
  sub <- coded.data[coded.data$Coder==Coders[i],]
  sub <- sub[!duplicated(sub$Document.ID,fromLast=T),]
  sub
})

data. <- do.call(rbind,data.)
coded.data <- data.

# ICR

X <- coded.data[,c("Document.ID","Coder", "RC.any","vic_blame", "empathy", "vic_trust","consent", "male.rape","critique" )]
colnames(X)[colnames(X) == "Document.ID"] <- 'AID'
X[is.na(X)] <- 0

# All
irr.mat.all <- icrFun(X=X,Coders=Coders,keepz=NULL,stripz=NULL,varz=NULL,nboot=1000,rnd=2)

